﻿using System;
using System.Collections.Generic;
using Allegro.MathInterface;

namespace Allegro.Mathlib
{
    public class ProjectedVectorFunction
    {

        public ProjectedVectorFunction(VectorFunction vf, int coordinateIndex)
        {
            SelectedCoordinateIndex = coordinateIndex;
            Func = vf;
        }
        /// <summary>
        /// Selected coordinate index
        /// </summary>
        public int SelectedCoordinateIndex { get; set; }

        /// <summary>
        /// Source vector function
        /// </summary>
        public VectorFunction Func { get; set; }

        /// <summary>
        /// When the function is evaluated, the vector function will be evaluated with argument set equal
        /// to vector, except for the selected coordinate, which will be assigned the value given by the function argument
        /// </summary>
        public double[] Vector { get; set; }

        /// <summary>
        /// Return the value of the vector function with argument equal to Vector, except for the selected coordinate, which 
        /// is assigned the value x.
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public double Value(double x)
        {
            Vector[SelectedCoordinateIndex] = x;
            return Func(Vector);
        }
    }
    public class Gradient : INumericalVectorDifferentiation
    {
        private ProjectedVectorFunction[] _projections;

        /// <summary>
        /// Construct a gradient using a specified differentiation method
        /// </summary>
        /// <param name="f"></param>
        /// <param name="diffMethod"></param>
        /// <param name="stepSizes"></param>
        public Gradient(VectorFunction f, IDifferentiation diffMethod, double[] stepSizes)
        {
            StepSizes = stepSizes;
            DifferentiationMethod = diffMethod;
            Func = f;
            int dim = stepSizes.Length;
            Dimension = dim;
            _projections = new ProjectedVectorFunction[dim];
            for (int i = 0; i < dim; i++)
            {
                _projections[i] = new ProjectedVectorFunction(f, i);
            }
        }
        /// <summary>
        /// Construct a gradient, using RiddersDifferentiationMethod as differentiation method
        /// </summary>
        /// <param name="f"></param>
        /// <param name="stepSizes"></param>
        public Gradient(VectorFunction f, double[] stepSizes)
            : this(f, new RiddersDifferentiationMethod(), stepSizes)
        {

        }
        public int Dimension { get; set; }


        #region Implementation of IVectorDifferentiation

        public double[] StepSizes { get; set; }
        public virtual IDifferentiation DifferentiationMethod { get; set; }

        public VectorFunction Func { get; set; }

        public void Df(double[] x, double[] dfdx)
        {
            int l = x.Length;

            for (int i = 0; i < l; i++)
            {
                // Store the original value of coordinate i.
                // For efficiency reasons, the differentiation method works directly on the input vector x, one coordinate
                // at a time.
                double x0 = x[i];

                _projections[i].Vector = x; // Assign the vector. The differentiation method will work on the input vector
                // Assign the function to be differentiated. This is the projection along coordinate i
                DifferentiationMethod.StepSize = StepSizes[i];
                DifferentiationMethod.Func = _projections[i].Value;
                // Differentiate
                dfdx[i] = DifferentiationMethod.Df(x0);

                // Restore the original value for coordinate i
                x[i] = x0;
            }
        }

        #endregion
    }
}
